% llrts_STAT.M
% compute posterior mean, sd, and confidence interval (HPD)
% based on Metropolis-Hastings output

% Taeyoung Doh
% created: 12/11/2007

%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for linear posterior draws
%	even integer (e.g., 2) for quadratic posterior draws
mhrun    = 0; 
mhrunid  = num2str(mhrun','%2.0f');

block1   = 1;		% starting block
blockn   = 10;		% ending block
nsim     = 10000;	% number of simulation draws in each block

hpdprob  = 0.90;	% probability for highest posterior density



%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;

% parameter names
pnames = feval(fnames_.pnames,msel);

% reset diary file
diaryfile = [path_.result 'stat_mhrun' mhrunid '.log'];
if exist(diaryfile,'file')
	delete(diaryfile);
end
diary(diaryfile);
diary off;

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------


% stind=1
% load path
 lpath = [path_.para 'mhrun' mhrunid  '/' ];
%  lpath = [path_.para 'mhrun' mhrunid '/' ];  
% initialize output
ndraws     = 0;
drawmean   = 0;
drawsqmean = 0;
sumdrawsq  = 0;
mpos       = zeros(blockn-block1+1,1); 
% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);

	% use only p.d. covariance matrix results for prior draws
	if mhrun==0
		parasim = parasim(find(retcode==0),:);
	end
     
	% number of simulations in each block
	nsimul = size(parasim,1);

	% collect simulations
	drawblock = parasim;
    drawdim   = size(drawblock,2);

	% compute sums of x(i) and x(i)^2 to be used to calculate means and s.d.
	drawmean   = drawmean + sum(drawblock,1);
	drawsqmean = drawsqmean + sum(drawblock.^2,1);
	sumdrawsq  = sumdrawsq + drawblock'*drawblock;
	ndraws     = ndraws + nsimul;
end

% compute means and standard deviations
drawmean   = drawmean/ndraws;
drawsqmean = drawsqmean/ndraws;
drawstdd   = sqrt(drawsqmean - drawmean.^2);
drawsig    = sumdrawsq/ndraws - drawmean'*drawmean;

%--------------------------------------------------------------------
% HPD Interval
%--------------------------------------------------------------------
drawci = zeros(2,drawdim);

for j=1:drawdim

	drawcol = [];
	ndraws  = 0;

	% loop for block
	for jblock=block1:blockn

		% load file
		lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
		load(lfile);

		% number of simulations in each block
		[nsimul,npara] = size(parasim);
        drawblock = parasim;  
        
		% Read only the j'th column
		drawcol = [drawcol ; drawblock(:,j)];
		ndraws  = ndraws + nsimul;

	end

	% calculate HPD for j'th column
	drawci(:,j) = hpdint(drawcol,hpdprob);

end


% save output
sfile = [path_.para 'stat_mhrun' mhrunid '.mat'];
save(sfile,'drawmean','drawstdd','drawci');

% print output 
fmt='%15.8f ';
disp(' ');
disp(' ');
disp(' ');
diary on;
disp(' ');
disp(sprintf('-------------------------------------------------------------------------'));
disp(         'Parameters      Mean            Stdev                   90% CI ');
for i=1:drawdim 
	disp(sprintf(['%s '  repmat(fmt,1,4) ],pnames(i,:),drawmean(i),drawstdd(i),drawci(2,i),drawci(1,i)));
end



